graficador <- function(a, b, l, f) {
    x <- seq(a, b, length.out = l)
    y <- x
    z <- outer(x, y, f)
    z[z < 0] = 0
    op <- par(bg = "white")
    persp(x, y, z, theta = 30, phi = 30, expand = 0.5,
        col = "lightblue", ltheta = 120, shade = 0.75,
        zlab = f, xlab = formalArgs(args(f))[1],
        ylab = formalArgs(args(f))[2]) -> res
}
tabla <- function(dataframe) 
  # Genera tablas con código de colores para un dataframe
  {
    gt(dataframe, rownames_to_stub = TRUE) %>%
        gt_color_rows(columns = dimnames.data.frame(dataframe)[[2]],
            palette = "ggsci::blue_material", domain = range(dataframe)) %>%
        gt_theme_nytimes()
}

FUNCION DE PRECIO función relu

price_function <- function(p_1, p_2, k = 400)
{
    return(1 - (exp(p_1/k)/(exp(p_1/k) + exp(p_2/k))))
}
graficador(3, 1000, 60, "price_function")

FUNCION DE DECISON

dicotomic_function_ps <- function(c_ps, c_xb, l_ps = 0.8) {
  return((1 - c_ps) * (1 - c_xb) * 0.75 +
           c_ps * (1 - c_xb) * 0.9 +
           c_xb * (1 - c_ps) * 0.3 +
           c_ps * c_xb * l_ps)
}

dicotomic_function_xb <- function(c_ps, c_xb, l_xb = 0.9) {
  return((1 - c_ps) * (1 - c_xb) * 0.7 +
           c_ps * (1 - c_xb) * 0.4 +
           c_xb * (1 - c_ps) * 0.8 +
           c_ps * c_xb * l_xb)
}

CONTROL DE PRECIO

price_control_ps <- function(p, m = 2000, c_ps, c_xb,
    l_ps = 0.8) {
    m = m * dicotomic_function_ps(c_ps, c_xb, l_ps = l_ps)
    return(0 + (p <= m) * (p>=0) *((-(p - m)/m)^2))
}
price_control_xb <- function(p, m = 2000, c_ps, c_xb,
    l_xb = 0.9) {
    m = m * dicotomic_function_xb(c_ps, c_xb, l_xb = l_xb)
    return(0 + (p <= m) * (p>=0) * (((p - m)/m)^2))
}

FUNCIONES DE VENTAS

ventas_ps <- function(p_ps, p_xb, c_ps = 1, c_xb = 0,
    k = 400, m = 2000, l_ps = 0.8)
    {
    return(5e+07 * dicotomic_function_ps(c_ps, c_xb,
        l_ps = l_ps) * price_function(p_ps, p_xb, k) *
        price_control_ps(p_ps, m, c_ps, c_xb, l_ps = l_ps))
}
ventas_xb <- function(p_ps, p_xb, c_ps = 0, c_xb = 0,
    k = 400, m = 2000, l_xb = 0.9)
    {
    return(5e+07 * dicotomic_function_xb(c_ps, c_xb,
        l_xb = l_xb) * price_function(p_xb, p_ps, k) *
        price_control_xb(p_xb, m, c_ps, c_xb, l_xb = l_xb))
}

FUNCIONES DE BENEFICIOS

beneficios_ps <- function(p_ps, p_xb, c_ps = 0, c_xb = 0,
    k = 400, m = 2000, l_ps = 0.8)
    {
    ventas_ps(p_ps, p_xb, c_ps = c_ps, c_xb = c_xb, k = k,
        m = m, l_ps = l_ps) * (p_ps - 100)
}

beneficios_xb <- function(p_ps, p_xb, c_ps = 0, c_xb = 0,
    k = 400, m = 2000, l_xb = 0.9)
    {
    ventas_xb(p_ps, p_xb, c_ps = c_ps, c_xb = c_xb, k = k,
        m = m, l_xb = l_xb) * (p_xb - 100)
}

graficador(1, 1000, 60, "beneficios_ps")

graficador(1, 1000, 60, "beneficios_xb")

FUNCIONES DE RESOLUCIÓN

# la siguiente función calcula el equilibrio de Nash
# del problema de Bertrand
resolucion  <- function(c_ps = 0,
                       c_xb = 0,
                       k = 400,
                       m = 2000,
                       x0 = 0,
                       y0 = 0,
                       N = 100,
                       l_ps = 0.8,
                       l_xb = 0.9,
                       eps=10^-4)
  # x0 e y0 corresponde al inicio
  # N corresponde al número de iteraciones
  # usaremos por defecto N = 100 por funcionar bien
{
  # generamos una matrix que contenga todas las iteraciones
  z = matrix(0, nrow = 2, ncol = N + 1)
  z[, 1] <- c(x0, y0)# almacenamos el punto inicial
  i = 0
  l=TRUE#criterio de parada por diverger
  aux2=TRUE #criterio de parada al converger
  # iniciamos el bucle
  while (all(i < N, aux2,l))
  {
    i = i + 1
    #almacenamos el óptimo para p_p
    aux = optim(z[1, i], function(x)
    {
      -beneficios_ps(x, z[2, i], c_ps,
                     c_xb, k, m, l_ps = l_ps)
    })$par
    #comprobamos que no sea demasiado alto
    if (aux <= (m * dicotomic_function_ps(c_ps, c_xb, l_ps = l_ps))) {
      z[1, i + 1] <- aux
    } else{
      l=FALSE
    }
    
    #almacenamos el óptimo para p_x
    
    aux = optim(z[2, i], function(y)
    {
      -beneficios_xb(z[1, i + 1], y,
                     c_ps, c_xb, k, m, l_xb = l_xb)
    })$par
    if (aux <= (m * dicotomic_function_xb(c_ps, c_xb, l_xb = l_xb))) {
      z[2, i + 1] <- aux
    } else{
      l=FALSE
    }
    aux2=dist(z[1, i + 1]-z[1, i])>eps
  }

  #nos devuelve los beneficios de cada empresa
  # y los precios de las consolas del equilibrio de Nash
  return(list(
    'beneficios' = c(
      beneficios_ps(
        p_ps = z[1,
                 i + 1],
        p_xb = z[2, i + 1],
        c_ps = c_ps,
        c_xb = c_xb,
        k,
        m,
        l_ps = l_ps
      ),
      beneficios_xb(
        p_ps = z[1, i + 1],
        p_xb = z[2, i + 1],
        c_ps = c_ps,
        c_xb = c_xb,
        k,
        m,
        l_xb = l_xb
      )
    ),
    'precios' = z[, i + 1],
    'converge' = l
  ))
}

SOLUCIÓN

# La siguiente función calcula la p y la q 
# de un problema 2x2
solucion2 <- function(e1,e2) {
# importamos las funciones de python
np <- import("numpy", convert = FALSE)
nash <- import("nashpy", convert = FALSE)
# creamos el juego
np1<-np$array(e1)

np2<-np$array(e2)

jg<- nash$Game(np1,np2)
#obtenemos el equilibrio 
equilibria <- jg$support_enumeration(non_degenerate = FALSE)
for (i in iterate(equilibria)){

}
#el último i contiene la estrategia mixta, si la hay
  return(list('p' = py_to_r(i)[[1]][1],
              'q' = py_to_r(i)[[2]][1]))
}
# esta función resuelve el problema completo
# devolviendonos la solución del mismo 

solucion <- function(k = 400,
                     m = 2000,
                     l_ps = 0.8,
                     l_xb = 0.9) {
  # calculamos los diferentes ejercicios de Bertrand
  precios1=matrix(
        nrow = 2,
        ncol = 2,
        dimnames = list(
          c("No hace crossplay", "Hace crossplay"),
          c("No hace crossplay", "Hace crossplay"))
        )
  precios2=precios1
  beneficios1=precios1
  beneficios2=precios1
  for (i in 0:1) {
    for (j in 0:1){
      r= resolucion(
        k = k,
        m = m,
        c_ps = j ,
        c_xb = i,
        l_ps = l_ps,
        l_xb = l_xb
      )
    precios1[j+1,i+1]=r$precios[1]
    precios2[j+1,i+1]=r$precios[2]
    beneficios1[j+1,i+1]=round(r$beneficios[1])
    beneficios2[j+1,i+1]=round(r$beneficios[2])
    }}
  # usamos el programa anterior para calcular 
  # la solución del problema 2x2
  solution2x2 = solucion2(precios1, precios2)
  p = solution2x2$p
  q = solution2x2$q
  # nos devuelve la solución junto con los 
  # beneficios esperados
  return(
    list(
      'p' = p,
      'q' = q,
      'beneficiosesperados' = list(
        'ps' = cbind(p, 1 - p) %*%
          beneficios1 %*%
          c(q, 1 - q),
        'xb' = cbind(p, 1 - p) %*%
          beneficios2 %*%
          c(q, 1 - q)
      ),
      'beneficiosps' = as.data.frame(beneficios1),
      'beneficiosxb' = as.data.frame(beneficios2),
      #obtenemos los precios que pone cada jugador en cada opción
      'preciosps' = as.data.frame(precios1),
      'preciosxb' = as.data.frame(precios2)
    )
  )
}
sol =solucion()
sol$beneficiosesperados
## $ps
##            [,1]
## [1,] 3665411948
## 
## $xb
##            [,1]
## [1,] 4376168454
sol$p
## [1] 0
sol$q
## [1] 0
sol$preciosps
##                   No hace crossplay Hace crossplay
## No hace crossplay          420.3572       249.2645
## Hace crossplay             440.2721       439.8871
sol$preciosxb
##                   No hace crossplay Hace crossplay
## No hace crossplay          408.0221       413.5474
## Hace crossplay             299.0318       460.4085
sol$beneficiosps
##                   No hace crossplay Hace crossplay
## No hace crossplay        3063835537      460007800
## Hace crossplay           3605449303     3665411948
sol$beneficiosxb
##                   No hace crossplay Hace crossplay
## No hace crossplay        2747973569     2749898042
## Hace crossplay            916862215     4376168454

ANÁLISIS DE SENSIBILIDAD

# Función que realiza un mallado 
# de la k y la m 
sensitivity <- function(lstk, lstm)
{
    #usamos las variables auxiliares 
    psolaux = c()# almacena las p obtenidas
    qsolaux = c()# almacena las q obtenidas
    pssolaux = c()# almacena los beneficios de 
    # ps obtenidos
    xbsolaux = c()# almacena los beneficios de 
    # xb obtenidos
    ps2solaux = c()# almacena los precios de 
    # ps obtenidos para c_ps=0 y c_xb=1
    xb2solaux = c()# almacena los precios de 
    # xb obtenidos para c_ps=0 y c_xb=1
    for (i in lstm)# bucle para la m
    {
        for (j in lstk)# bucle para la k
        {   # obtenemos la solución
            sol = solucion(k = j, m = i)
            # almacenamos los datos que interesan
            psolaux = c(psolaux, round(sol$p,
                4))
            qsolaux = c(qsolaux, round(sol$q,
                4))
            pssolaux = c(pssolaux,
                round(sol$beneficiosesperados$ps,
                  4))
            xbsolaux = c(xbsolaux,
                round(sol$beneficiosesperados$xb,
                  4))
            ps2solaux = c(ps2solaux,
                round(sol$preciosps[2,
                  2], 2))
            xb2solaux = c(xb2solaux,
                round(sol$preciosxb[2,
                  2], 2))
        }
    }
    # transformamos los datos a dataframes
    plista = as.data.frame(matrix(psolaux,
        nrow = length(lstk), ncol = length(lstm),
        dimnames = list(paste("k=",
            lstk), paste("m=", lstm))))
    qlista = as.data.frame(matrix(qsolaux,
        nrow = length(lstk), ncol = length(lstm),
        dimnames = list(paste("k=",
            lstk), paste("m=", lstm))))
    pslista = as.data.frame(matrix(pssolaux,
        nrow = length(lstk), ncol = length(lstm),
        dimnames = list(paste("k=",
            lstk), paste("m=", lstm))))
    xblista = as.data.frame(matrix(xbsolaux,
        nrow = length(lstk), ncol = length(lstm),
        dimnames = list(paste("k=",
            lstk), paste("m=", lstm))))
    xb2lista = as.data.frame(matrix(xb2solaux,
        nrow = length(lstk), ncol = length(lstm),
        dimnames = list(paste("k=",
            lstk), paste("m=", lstm))))
    ps2lista = as.data.frame(matrix(ps2solaux,
        nrow = length(lstk), ncol = length(lstm),
        dimnames = list(paste("k=",
            lstk), paste("m=", lstm))))
    # devolvemos los diferentes dataframes
    return(list('p' = plista, 'q' = qlista,
        'beneficiosps' = pslista, 'beneficiosxb' = xblista,
        'preciosps' = ps2lista, 'preciosxb' = xb2lista))
}
# creamos el mallado
lstk = seq(200, 600, length.out = 9)
lstm = seq(1000, 3000, length.out = 9)
# evaluamos el mallado 
lista = sensitivity(lstk, lstm)
# obtenemos todas las tablas

tabla(lista[["p"]])
m= 1000 m= 1250 m= 1500 m= 1750 m= 2000 m= 2250 m= 2500 m= 2750 m= 3000
k= 200 0 0 0 0 0.8752 0.848 0.8250 0.8056 0.7898
k= 250 0 0 0 0 0.0000 0.000 0.8726 0.8519 0.8328
k= 300 0 0 0 0 0.0000 0.000 0.0000 0.0000 0.8721
k= 350 0 0 0 0 0.0000 0.000 0.0000 0.0000 0.0000
k= 400 0 0 0 0 0.0000 0.000 0.0000 0.0000 0.0000
k= 450 0 0 0 0 0.0000 0.000 0.0000 0.0000 0.0000
k= 500 0 0 0 0 0.0000 0.000 0.0000 0.0000 0.0000
k= 550 0 0 0 0 0.0000 0.000 0.0000 0.0000 0.0000
k= 600 0 0 0 0 0.0000 0.000 0.0000 0.0000 0.0000
tabla(lista[["q"]])
m= 1000 m= 1250 m= 1500 m= 1750 m= 2000 m= 2250 m= 2500 m= 2750 m= 3000
k= 200 0 0 0 0 0.9928 0.977 0.9651 0.9548 0.9470
k= 250 0 0 0 0 0.0000 0.000 0.9920 0.9796 0.9689
k= 300 0 0 0 0 0.0000 0.000 0.0000 0.0000 0.9917
k= 350 0 0 0 0 0.0000 0.000 0.0000 0.0000 0.0000
k= 400 0 0 0 0 0.0000 0.000 0.0000 0.0000 0.0000
k= 450 0 0 0 0 0.0000 0.000 0.0000 0.0000 0.0000
k= 500 0 0 0 0 0.0000 0.000 0.0000 0.0000 0.0000
k= 550 0 0 0 0 0.0000 0.000 0.0000 0.0000 0.0000
k= 600 0 0 0 0 0.0000 0.000 0.0000 0.0000 0.0000
tabla(lista[['beneficiosps']])
m= 1000 m= 1250 m= 1500 m= 1750 m= 2000 m= 2250 m= 2500 m= 2750 m= 3000
k= 200 1507584083 1965591008 2380143772 2752750928 2652873864 2906595062 3140921257 3354637132 3551799715
k= 250 1546131330 2041127306 2501788127 2926159824 3314630443 3669686131 3444385258 3695436686 3929660465
k= 300 1568107579 2087385417 2579991822 3041689354 3472337098 3871936572 4241712304 4585411821 4236637708
k= 350 1581317042 2116865447 2631766986 3120899396 3583849435 4018537734 4426469866 4808418765 5165654640
k= 400 1589610198 2136159879 2667660319 3178076864 3665411948 4127766690 4566005972 4979761959 5370182729
k= 450 1594821540 2149494742 2692764130 3218917373 3725420728 4210351901 4672534439 5113045012 5532264982
k= 500 1598154405 2158630791 2710735907 3249232976 3770762571 4273313057 4755872762 5218676287 5660733657
k= 550 1600444871 2164962963 2724059154 3271892288 3805303363 4322349683 4820988533 5302190893 5765104296
k= 600 1601867661 2169741342 2733853579 3289212343 3832239148 4360797813 4874107560 5369454961 5849185658
tabla(lista[['beneficiosxb']])
m= 1000 m= 1250 m= 1500 m= 1750 m= 2000 m= 2250 m= 2500 m= 2750 m= 3000
k= 200 1839819353 2361427842 2830308934 3249545230 2194850974 2395607578 2577488075 2744601586 2898366183
k= 250 1901269854 2468318160 2991746669 3472061795 3909612095 4308629825 2858576629 3059381206 3243179154
k= 300 1941150950 2539320374 3101782734 3626632952 4112964205 4564244623 4979848153 5365708599 3527160129
k= 350 1968901134 2588896408 3179063504 3737922487 4262993529 4754425596 5215309730 5645356619 6046917677
k= 400 1989055132 2625235912 3237404138 3821430326 4376168454 4900972430 5397094732 5864128281 6304138902
k= 450 2004316972 2652984610 3281296274 3885561098 4464085862 5015532028 5540123887 6039106288 6512146528
k= 500 2016297824 2674562585 3315837431 3936115123 4533714984 5106978440 5655492205 6180745235 6681047097
k= 550 2025885318 2691853692 3343519388 3976700336 4589626178 5181138827 5749669826 6296532312 6821040583
k= 600 2033627763 2706077134 3366207940 4010024095 4635529465 5242055425 5827721737 6392924016 6938014217
tabla(lista[['preciosps']])
m= 1000 m= 1250 m= 1500 m= 1750 m= 2000 m= 2250 m= 2500 m= 2750 m= 3000
k= 200 262.54 290.32 312.87 331.48 346.95 359.97 371.05 380.59 388.87
k= 250 274.02 306.94 334.62 358.20 378.23 395.51 410.39 423.46 434.87
k= 300 282.32 319.34 351.29 379.04 403.09 424.29 442.75 459.27 473.89
k= 350 288.62 328.83 364.15 395.57 423.31 447.84 469.87 489.49 507.09
k= 400 293.51 336.33 374.72 409.05 439.89 467.57 492.61 515.12 535.52
k= 450 297.44 342.45 383.21 420.17 453.75 484.18 511.85 537.17 560.23
k= 500 300.68 347.47 390.30 429.49 465.45 498.32 528.41 556.22 581.63
k= 550 303.40 351.69 396.29 437.39 475.37 510.47 542.75 572.72 600.44
k= 600 305.64 355.32 401.43 444.19 483.92 520.97 555.28 587.21 617.02
tabla(lista[['preciosxb']])
m= 1000 m= 1250 m= 1500 m= 1750 m= 2000 m= 2250 m= 2500 m= 2750 m= 3000
k= 200 273.54 300.44 322.04 339.78 354.46 366.76 377.25 386.25 394.13
k= 250 287.86 320.31 347.22 369.85 389.01 405.43 419.60 431.91 442.73
k= 300 298.52 335.58 367.05 393.95 417.26 437.50 455.14 470.86 484.71
k= 350 306.74 347.59 382.89 413.54 440.64 464.39 485.49 504.27 521.04
k= 400 313.27 357.20 395.95 430.06 460.41 487.34 511.54 533.22 552.76
k= 450 318.52 365.17 406.73 443.80 477.12 507.07 533.99 558.45 580.76
k= 500 322.81 371.77 415.78 455.54 491.49 524.07 553.66 580.73 605.37
k= 550 326.50 377.29 423.56 465.59 503.91 538.93 570.82 600.27 627.40
k= 600 329.59 382.16 430.22 474.35 514.80 551.95 586.29 617.64 646.92
# Función que realiza un mallado
# de la l_ps y la l_xb
sensitivity2 <- function(n,a,b)
{
  lstps=rbeta(n,a,b)
      #usamos las variables auxiliares
    psolaux = c()# almacena las p obtenidas
    qsolaux = c()# almacena las q obtenidas
    pssolaux = c()# almacena los beneficios de
    # ps obtenidos
    xbsolaux = c()# almacena los beneficios de
    # xb obtenidos
    ps2solaux = c()# almacena los precios de
    # ps obtenidos para c_ps=0 y c_xb=1
    xb2solaux = c()# almacena los precios de
    # xb obtenidos para c_ps=0 y c_xb=1
    for (l_ps in lstps)# bucle para la l_ps
        {   # obtenemos la solución
            sol = solucion(l_ps = l_ps)
            # almacenamos los datos que interesan
            psolaux = c(psolaux, round(sol$p,
                4))
            qsolaux = c(qsolaux, round(sol$q,
                4))
            pssolaux = c(pssolaux,
                round(sol$beneficiosesperados$ps,
                  4))
            xbsolaux = c(xbsolaux,
                round(sol$beneficiosesperados$xb,
                  4))
            ps2solaux = c(ps2solaux,
                round(sol$preciosps[2,
                  2], 2))
            xb2solaux = c(xb2solaux,
                round(sol$preciosxb[2,
                  2], 2))
        }

    
    # devolvemos los diferentes dataframes
    return(list('p' = psolaux, 'q' = qsolaux,
        'beneficiosps' = pssolaux, 'beneficiosxb' = xbsolaux,
        'preciosps' = ps2solaux, 'preciosxb' = xb2solaux))
}
set.seed(2)
n=100
a=8
b=2
lista2 = sensitivity2(n,a,b)
#obtenemos todas las tablas

mean(lista2[["p"]])
## [1] 0
mean(lista2[["q"]])
## [1] 0
histo<-function(lst){
  m=mean(lst)
  mo=mlv(lst,method = 'shorth')
  
  q1=quantile(lst, 0.025)
  q2=quantile(lst, 0.975)
  
  hist(lst,freq=FALSE)
  abline(v =m,col = "blue",lwd = 2)
  abline(v =q1,col = "red",lwd = 2)
  abline(v =q2,col = "red",lwd = 2)
  abline(v =mo,col = "green",lwd = 2)
  return(paste(
    'El cuantil para 0.025 es ' , as.character(q1),'
    El cuantil para 0.975 es ',as.character(q2),'
    La media es ' ,as.character(m),'
    La moda es ' ,as.character(mo),'
    La desviación típica es ',as.character(sd(lst))
  ))
}
# analizamos los beneficios de Playstation
histo(lista2[['beneficiosps']])

## [1] "El cuantil para 0.025 es  1615273846.525 \n    El cuantil para 0.975 es  5198874991.925 \n    La media es  3802338582.66 \n    La moda es  4302440780.54 \n    La desviación típica es  1006283136.00399"
# analizamos los beneficios de Xbox
histo(lista2[['beneficiosxb']])

## [1] "El cuantil para 0.025 es  3883555405.2 \n    El cuantil para 0.975 es  4605210455.425 \n    La media es  4374980663.61 \n    La moda es  4482525143.96 \n    La desviación típica es  206616743.334366"
histo(lista2[['preciosps']])

## [1] "El cuantil para 0.025 es  350.83675 \n    El cuantil para 0.975 es  480.531 \n    La media es  439.4743 \n    La moda es  458.7956 \n    La desviación típica es  37.2257921408686"
histo(lista2[['preciosxb']])

## [1] "El cuantil para 0.025 es  447.37075 \n    El cuantil para 0.975 es  466.60975 \n    La media es  460.3759 \n    La moda es  463.4226 \n    La desviación típica es  5.48209631344607"